home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / sspmv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  8.1 KB  |  266 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE SSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
  5. *     .. Scalar Arguments ..
  6.       REAL               ALPHA, BETA
  7.       INTEGER            INCX, INCY, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       REAL               AP( * ), X( * ), Y( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  SSPMV  performs the matrix-vector operation
  17. *
  18. *     y := alpha*A*x + beta*y,
  19. *
  20. *  where alpha and beta are scalars, x and y are n element vectors and
  21. *  A is an n by n symmetric matrix, supplied in packed form.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the matrix A is supplied in the packed
  29. *           array AP as follows:
  30. *
  31. *              UPLO = 'U' or 'u'   The upper triangular part of A is
  32. *                                  supplied in AP.
  33. *
  34. *              UPLO = 'L' or 'l'   The lower triangular part of A is
  35. *                                  supplied in AP.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - REAL            .
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  AP     - REAL             array of DIMENSION at least
  49. *           ( ( n*( n + 1 ) )/2 ).
  50. *           Before entry with UPLO = 'U' or 'u', the array AP must
  51. *           contain the upper triangular part of the symmetric matrix
  52. *           packed sequentially, column by column, so that AP( 1 )
  53. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  54. *           and a( 2, 2 ) respectively, and so on.
  55. *           Before entry with UPLO = 'L' or 'l', the array AP must
  56. *           contain the lower triangular part of the symmetric matrix
  57. *           packed sequentially, column by column, so that AP( 1 )
  58. *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  59. *           and a( 3, 1 ) respectively, and so on.
  60. *           Unchanged on exit.
  61. *
  62. *  X      - REAL             array of dimension at least
  63. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  64. *           Before entry, the incremented array X must contain the n
  65. *           element vector x.
  66. *           Unchanged on exit.
  67. *
  68. *  INCX   - INTEGER.
  69. *           On entry, INCX specifies the increment for the elements of
  70. *           X. INCX must not be zero.
  71. *           Unchanged on exit.
  72. *
  73. *  BETA   - REAL            .
  74. *           On entry, BETA specifies the scalar beta. When BETA is
  75. *           supplied as zero then Y need not be set on input.
  76. *           Unchanged on exit.
  77. *
  78. *  Y      - REAL             array of dimension at least
  79. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  80. *           Before entry, the incremented array Y must contain the n
  81. *           element vector y. On exit, Y is overwritten by the updated
  82. *           vector y.
  83. *
  84. *  INCY   - INTEGER.
  85. *           On entry, INCY specifies the increment for the elements of
  86. *           Y. INCY must not be zero.
  87. *           Unchanged on exit.
  88. *
  89. *
  90. *  Level 2 Blas routine.
  91. *
  92. *  -- Written on 22-October-1986.
  93. *     Jack Dongarra, Argonne National Lab.
  94. *     Jeremy Du Croz, Nag Central Office.
  95. *     Sven Hammarling, Nag Central Office.
  96. *     Richard Hanson, Sandia National Labs.
  97. *
  98. *
  99. *     .. Parameters ..
  100.       REAL               ONE         , ZERO
  101.       PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
  102. *     .. Local Scalars ..
  103.       REAL               TEMP1, TEMP2
  104.       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY
  105. *     .. External Functions ..
  106.       LOGICAL            LSAME
  107.       EXTERNAL           LSAME
  108. *     .. External Subroutines ..
  109.       EXTERNAL           XERBLA
  110. *     ..
  111. *     .. Executable Statements ..
  112. *
  113. *     Test the input parameters.
  114. *
  115.       INFO = 0
  116.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  117.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  118.          INFO = 1
  119.       ELSE IF( N.LT.0 )THEN
  120.          INFO = 2
  121.       ELSE IF( INCX.EQ.0 )THEN
  122.          INFO = 6
  123.       ELSE IF( INCY.EQ.0 )THEN
  124.          INFO = 9
  125.       END IF
  126.       IF( INFO.NE.0 )THEN
  127.          CALL XERBLA( 'SSPMV ', INFO )
  128.          RETURN
  129.       END IF
  130. *
  131. *     Quick return if possible.
  132. *
  133.       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  134.      $   RETURN
  135. *
  136. *     Set up the start points in  X  and  Y.
  137. *
  138.       IF( INCX.GT.0 )THEN
  139.          KX = 1
  140.       ELSE
  141.          KX = 1 - ( N - 1 )*INCX
  142.       END IF
  143.       IF( INCY.GT.0 )THEN
  144.          KY = 1
  145.       ELSE
  146.          KY = 1 - ( N - 1 )*INCY
  147.       END IF
  148. *
  149. *     Start the operations. In this version the elements of the array AP
  150. *     are accessed sequentially with one pass through AP.
  151. *
  152. *     First form  y := beta*y.
  153. *
  154.       IF( BETA.NE.ONE )THEN
  155.          IF( INCY.EQ.1 )THEN
  156.             IF( BETA.EQ.ZERO )THEN
  157.                DO 10, I = 1, N
  158.                   Y( I ) = ZERO
  159.    10          CONTINUE
  160.             ELSE
  161.                DO 20, I = 1, N
  162.                   Y( I ) = BETA*Y( I )
  163.    20          CONTINUE
  164.             END IF
  165.          ELSE
  166.             IY = KY
  167.             IF( BETA.EQ.ZERO )THEN
  168.                DO 30, I = 1, N
  169.                   Y( IY ) = ZERO
  170.                   IY      = IY   + INCY
  171.    30          CONTINUE
  172.             ELSE
  173.                DO 40, I = 1, N
  174.                   Y( IY ) = BETA*Y( IY )
  175.                   IY      = IY           + INCY
  176.    40          CONTINUE
  177.             END IF
  178.          END IF
  179.       END IF
  180.       IF( ALPHA.EQ.ZERO )
  181.      $   RETURN
  182.       KK = 1
  183.       IF( LSAME( UPLO, 'U' ) )THEN
  184. *
  185. *        Form  y  when AP contains the upper triangle.
  186. *
  187.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  188.             DO 60, J = 1, N
  189.                TEMP1 = ALPHA*X( J )
  190.                TEMP2 = ZERO
  191.                K     = KK
  192.                DO 50, I = 1, J - 1
  193.                   Y( I ) = Y( I ) + TEMP1*AP( K )
  194.                   TEMP2  = TEMP2  + AP( K )*X( I )
  195.                   K      = K      + 1
  196.    50          CONTINUE
  197.                Y( J ) = Y( J ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2
  198.                KK     = KK     + J
  199.    60       CONTINUE
  200.          ELSE
  201.             JX = KX
  202.             JY = KY
  203.             DO 80, J = 1, N
  204.                TEMP1 = ALPHA*X( JX )
  205.                TEMP2 = ZERO
  206.                IX    = KX
  207.                IY    = KY
  208.                DO 70, K = KK, KK + J - 2
  209.                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
  210.                   TEMP2   = TEMP2   + AP( K )*X( IX )
  211.                   IX      = IX      + INCX
  212.                   IY      = IY      + INCY
  213.    70          CONTINUE
  214.                Y( JY ) = Y( JY ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2
  215.                JX      = JX      + INCX
  216.                JY      = JY      + INCY
  217.                KK      = KK      + J
  218.    80       CONTINUE
  219.          END IF
  220.       ELSE
  221. *
  222. *        Form  y  when AP contains the lower triangle.
  223. *
  224.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  225.             DO 100, J = 1, N
  226.                TEMP1  = ALPHA*X( J )
  227.                TEMP2  = ZERO
  228.                Y( J ) = Y( J )       + TEMP1*AP( KK )
  229.                K      = KK           + 1
  230.                DO 90, I = J + 1, N
  231.                   Y( I ) = Y( I ) + TEMP1*AP( K )
  232.                   TEMP2  = TEMP2  + AP( K )*X( I )
  233.                   K      = K      + 1
  234.    90          CONTINUE
  235.                Y( J ) = Y( J ) + ALPHA*TEMP2
  236.                KK     = KK     + ( N - J + 1 )
  237.   100       CONTINUE
  238.          ELSE
  239.             JX = KX
  240.             JY = KY
  241.             DO 120, J = 1, N
  242.                TEMP1   = ALPHA*X( JX )
  243.                TEMP2   = ZERO
  244.                Y( JY ) = Y( JY )       + TEMP1*AP( KK )
  245.                IX      = JX
  246.                IY      = JY
  247.                DO 110, K = KK + 1, KK + N - J
  248.                   IX      = IX      + INCX
  249.                   IY      = IY      + INCY
  250.                   Y( IY ) = Y( IY ) + TEMP1*AP( K )
  251.                   TEMP2   = TEMP2   + AP( K )*X( IX )
  252.   110          CONTINUE
  253.                Y( JY ) = Y( JY ) + ALPHA*TEMP2
  254.                JX      = JX      + INCX
  255.                JY      = JY      + INCY
  256.                KK      = KK      + ( N - J + 1 )
  257.   120       CONTINUE
  258.          END IF
  259.       END IF
  260. *
  261.       RETURN
  262. *
  263. *     End of SSPMV .
  264. *
  265.       END
  266.